 # coding=gb18030


import xlrd
import numpy
import openpyxl
import time
from gurobipy import *

g_node_list = []
g_link_list = []
g_agent_list = []
g_state_list = []


g_number_of_nodes = 0
g_number_of_links = 0
g_number_of_agents = 0
g_number_of_states = 0

# need to pre determined
g_number_of_time_intervals = 120 + 1
g_number_of_od_pairs = 2 + 1

     
class Node:
    def __init__(self):
        self.node_id = 0
        self.type = 0
        self.ingoing_node_list = []
        self.ingoing_link_list = []
        self.outgoing_node_list = []
        self.outgoing_link_list = []
        self.cumulative_departure_count = []
        self.cumulative_arrival_count = []


class Link:
    def __init__(self):
        self.link_id = 0
        self.from_node_id = 0
        self.to_node_id = 0
        self.type = 0
        self.cost = 0
        self.free_flow_link_travel_time = 0
        self.cost = 0
        self.link_capacity = 0
        self.node_sequence = []
        self.time_sequence= []
        self.flow_sequence = []
        self.link_sequence = []
        


class Agent:
    def __init__(self):
        self.agent_id = 0
        self.departure_time = 0
        self.arrival_time = 0
        self.from_node_id = 0
        self.to_node_id = 0
        self.reservation_flag = 0
        self.departure_time_for_non_reservation_agent = 0
        self.link_sequence_in_path = []
        self.node_sequence = []
        self.time_sequence = []
        self.link_sequence = []
        self.time_sequence = []
        self.node_sequence = []
        self.optimized_departure_time = 0
        self.optimized_arrival_time = 0

        
  
def g_read_input_data():    
    
    #initialization

    global g_number_of_nodes
    global g_number_of_links
    global g_number_of_agents


    # read node
    file = xlrd.open_workbook("input_node.xlsx") # open file
    sheet = file.sheet_by_index(0) # open sheet
    node = Node()
    node.node_id = 0
    g_node_list.append(node)
    g_number_of_nodes += 1    
    for row in range(1, sheet.nrows):
        try:
            node = Node()
            node.node_id = int(sheet.cell_value(row, 0))
            node.type = int(sheet.cell_value(row, 1))
            g_node_list.append(node)
            g_number_of_nodes += 1
            print('reading {} passenger nodes..'.format(g_number_of_nodes))
        except:
            print('Read error. Check your passenger node file')
    print('passenger nodes_number:{}'.format(g_number_of_nodes))  
    
    # read links
    file = xlrd.open_workbook("input_link.xlsx") # open file
    sheet = file.sheet_by_index(0) # open sheet
    link = Link()
    link.link_id = 0
    g_link_list.append(link)
    g_number_of_links += 1
    for row in range(1, sheet.nrows):
        try:
            link = Link()
            link.link_id = int(sheet.cell_value(row, 0))
            link.from_node_id = int(sheet.cell_value(row, 1))
            link.to_node_id = int(sheet.cell_value(row, 2))
            link.type = int(sheet.cell_value(row, 3))
            link.free_flow_link_travel_time = int(sheet.cell_value(row, 4))
            link.cost = float(sheet.cell_value(row, 5))
            link.link_capacity = int(sheet.cell_value(row, 6))
            g_link_list.append(link)
            g_number_of_links += 1
            print('reading {} passenger links..'.format(g_number_of_links))
        except:
            print('Read error. Check your passenger link file')
    print('passenger links_number:{}'.format(g_number_of_links))
    
    # read buses
    file = xlrd.open_workbook("input_agent.xlsx") # open file
    sheet = file.sheet_by_index(0) # open sheet
    agent = Agent()
    agent.agent_id = 0
    g_agent_list.append(agent)
    g_number_of_agents += 1
    for row in range(1, sheet.nrows):
        try:
            agent = Agent()
            agent.agent_id = int(sheet.cell_value(row, 0))
            agent.departure_time = int(sheet.cell_value(row, 1))
            agent.arrival_time = int(sheet.cell_value(row, 2))
            agent.from_node_id = int(sheet.cell_value(row, 3))
            agent.to_node_id = int(sheet.cell_value(row, 4))
            agent.reservation_flag = int(sheet.cell_value(row, 5))
            agent.departure_time_for_non_reservation_agent = int(sheet.cell_value(row, 6))
            link_sequence_in_path_1 = str(sheet.cell_value(row, 7))
            link_sequence_in_path_2 = link_sequence_in_path_1.strip().split(';')
            agent.link_sequence_in_path = [int (link) for link in link_sequence_in_path_2]
            g_agent_list.append(agent)
            g_number_of_agents += 1
            print('reading {} agents..'.format(g_number_of_agents))
        except:
            print('Read. Check your passenger file')
    print('passengers_number:{}'.format(g_number_of_agents))  
    
    

              
                         
def  g_generate_in_out_going_link():   

    # record ingoing and outgoing nodes and links for each passenger node       
    for l in range(1, g_number_of_links):
        link_id = g_link_list[l].link_id
        from_node_id = g_link_list[l].from_node_id
        to_node_id = g_link_list[l].to_node_id     
        g_node_list[to_node_id].ingoing_link_list.append(link_id)
        g_node_list[to_node_id].ingoing_node_list.append(from_node_id)
        g_node_list[from_node_id].outgoing_link_list.append(link_id)
        g_node_list[from_node_id].outgoing_node_list.append(to_node_id)
        
    return()


def g_use_Gurobi_to_sovle_the_model():
    
    global node_set
    global link_set
    global time_set
    global agent_set
    global x
    global y
    global z
    global u
    global departure_before
    global departure_after
    global arrival_before
    global arrival_after
    global non_reservation_arrival_before
    global non_reservation_arrival_after
    global FIFO_link_set
    
    
   
    departure_before = 0
    departure_after = 0
    arrival_before = 120
    arrival_after = 120
    
    non_reservation_arrival_before = 0
    non_reservation_arrival_after = 0

       
        
    m = Model("agent_reservation_model")
    
    # establish sets
    node_set = []
    link_set = []
    time_set = []
    agent_set = []
    
    for n in range(1, g_number_of_nodes):
        node_id = g_node_list[n].node_id
        node_set.append(node_id)    
            
    for l in range(1, g_number_of_links):
        link_id = g_link_list[l].link_id
        link_set.append(link_id)
           
    for t in range(1, g_number_of_time_intervals):
        time_set.append(t)
        
    for a in range(1, g_number_of_agents):
        agent_id = g_agent_list[a].agent_id
        agent_set.append(agent_id)
    
    # decision variables
    x = m.addVars(agent_set, link_set, time_set, time_set, vtype = GRB.BINARY, name = 'x')
    y = m.addVars(link_set, time_set, lb = 0, ub = 100, vtype = GRB.CONTINUOUS, name = 'y')
    z = m.addVars(agent_set, agent_set, link_set, vtype = GRB.BINARY, name = 'z')
    
    print("finish decision variable formulation")
    
  
    
    objective_function = sum(sum(sum(sum(x[agent_id, link_id, from_t, to_t] * g_link_list[link_id].free_flow_link_travel_time       \
                         for to_t in time_set if (to_t == from_t + g_link_list[link_id].free_flow_link_travel_time)) for from_t in time_set) for link_id in link_set) for agent_id in agent_set)
    print("finish objective function formulation")
    m.setObjective(objective_function, GRB.MINIMIZE)
    
    # infeasible solution space
    for agent_id in agent_set:
        for link_id in link_set:
            if (g_link_list[link_id].type == 1):
                for from_t in time_set:
                     to_t = from_t + g_link_list[link_id].free_flow_link_travel_time
                     for other_t in time_set:
                         if (other_t != to_t):
                             m.addConstr(x[agent_id, link_id, from_t, other_t] == 0, "c1")
                                          
    print("finish infeasible solution space formulation")
    
    # agent flow balance constraint
    # origin constraint
    for agent_id in agent_set:
        # passenger without reservation
        if (g_agent_list[agent_id].reservation_flag == 0):
            from_node_id = g_agent_list[agent_id].from_node_id
            to_node_id = g_agent_list[agent_id].to_node_id
            departure_time = g_agent_list[agent_id].departure_time
            arrival_time = g_agent_list[agent_id].arrival_time
            departure_time_for_non_reservation_agent = g_agent_list[agent_id].departure_time_for_non_reservation_agent
            m.addConstr(sum(sum(x[agent_id, outgoing_link_id, departure_time_for_non_reservation_agent, to_t] for to_t in time_set if (to_t == departure_time_for_non_reservation_agent + g_link_list[outgoing_link_id].free_flow_link_travel_time)) for outgoing_link_id in g_node_list[from_node_id].outgoing_link_list) == 1, "c2")
        if (g_agent_list[agent_id].reservation_flag == 1):
            from_node_id = g_agent_list[agent_id].from_node_id
            to_node_id = g_agent_list[agent_id].to_node_id
            departure_time = g_agent_list[agent_id].departure_time
            arrival_time = g_agent_list[agent_id].arrival_time
            departure_time_for_non_reservation_agent = g_agent_list[agent_id].departure_time_for_non_reservation_agent
            # establish departure time window
            departure_time_window = []
            for t in range(max(1, departure_time - departure_before), min(departure_time + departure_after + 1, g_number_of_time_intervals)):
            # for t in range(1, g_number_of_time_intervals):
                departure_time_window.append(t)
            m.addConstr(sum(sum(sum(x[agent_id, outgoing_link_id, scheduled_departure_time, to_t] for to_t in time_set if (to_t == scheduled_departure_time + g_link_list[outgoing_link_id].free_flow_link_travel_time)) for scheduled_departure_time in departure_time_window) for outgoing_link_id in g_node_list[from_node_id].outgoing_link_list) == 1, "c2")
    print("finish origin flow balance formulation")
    
    # destination constraint
    for agent_id in agent_set:
        if (g_agent_list[agent_id].reservation_flag == 0):
            from_node_id = g_agent_list[agent_id].from_node_id
            to_node_id = g_agent_list[agent_id].to_node_id
            departure_time = g_agent_list[agent_id].departure_time
            arrival_time = g_agent_list[agent_id].arrival_time
            # establish departure time window
            arrival_time_window = []
            for t in range(max(1, arrival_time - non_reservation_arrival_before), min(arrival_time + non_reservation_arrival_after + 1, g_number_of_time_intervals)):
            # for t in range(1, arrival_time + 1):
                arrival_time_window.append(t)
            m.addConstr(sum(sum(sum(x[agent_id, ingoing_link_id, from_t, arrival_time] for from_t in time_set if (from_t == arrival_time - g_link_list[ingoing_link_id].free_flow_link_travel_time)) for arrival_time in arrival_time_window) for ingoing_link_id in g_node_list[to_node_id].ingoing_link_list) == 1, "c3")        
        if (g_agent_list[agent_id].reservation_flag == 1):
            from_node_id = g_agent_list[agent_id].from_node_id
            to_node_id = g_agent_list[agent_id].to_node_id
            departure_time = g_agent_list[agent_id].departure_time
            arrival_time = g_agent_list[agent_id].arrival_time
            # establish departure time window
            arrival_time_window = []
            for t in range(max(1, arrival_time - arrival_before), min(arrival_time + arrival_after + 1, g_number_of_time_intervals)):
            # for t in range(1, arrival_time + 1):
                arrival_time_window.append(t)
            m.addConstr(sum(sum(sum(x[agent_id, ingoing_link_id, from_t, arrival_time] for from_t in time_set if (from_t == arrival_time - g_link_list[ingoing_link_id].free_flow_link_travel_time)) for arrival_time in arrival_time_window) for ingoing_link_id in g_node_list[to_node_id].ingoing_link_list) == 1, "c3")
    print("finish destination flow balance formulation")
    
    # intermidate vertex constraint
    for agent_id in agent_set:
        from_node_id = g_agent_list[agent_id].from_node_id
        to_node_id = g_agent_list[agent_id].to_node_id
        departure_time = g_agent_list[agent_id].departure_time
        arrival_time = g_agent_list[agent_id].arrival_time
        for inter_node_id in node_set:
            if (g_node_list[inter_node_id].type == 1):
                for inter_t in time_set:
                    m.addConstr(sum(sum(x[agent_id, ingoing_link_id, from_t, inter_t] for from_t in time_set if (from_t == inter_t - g_link_list[ingoing_link_id].free_flow_link_travel_time)) for ingoing_link_id in g_node_list[inter_node_id].ingoing_link_list) ==        \
                                sum(sum(x[agent_id, outgoing_link_id, inter_t, to_t] for to_t in time_set if (to_t == inter_t + g_link_list[outgoing_link_id].free_flow_link_travel_time)) for outgoing_link_id in g_node_list[inter_node_id].outgoing_link_list), "c4")
    print("finish intermidate flow balance formulation")

    # link capacity constraint
    for link_id in link_set:
        link_capacity = g_link_list[link_id].link_capacity
        for from_t in time_set:
            to_t = from_t + g_link_list[link_id].free_flow_link_travel_time
            if (to_t in time_set):
                m.addConstr(sum(x[agent_id, link_id, from_t, to_t] for agent_id in agent_set) <= link_capacity, "c5")
    print("finish link capacity constraint")
    
                    
    # calculate time_dependent link volume
    for link_id in link_set:
        for from_t in time_set:
            to_t = from_t + g_link_list[link_id].free_flow_link_travel_time
            if (to_t in time_set):
                m.addConstr(sum(x[agent_id, link_id, from_t, to_t] for agent_id in agent_set) == y[link_id, from_t], "c6")
    print("finish time_dependent link volume constraint")
    
    # first-in-first-out constraint
    for link_id in link_set:
        if (g_link_list[link_id].type == 1):
            for agent_id in agent_set:
                for another_agent_id in agent_set:
                    if (link_id in g_agent_list[agent_id].link_sequence_in_path) and (link_id in g_agent_list[another_agent_id].link_sequence_in_path):
                        m.addConstr(sum(sum(x[agent_id, link_id, from_t, to_t] * from_t for to_t in time_set if (to_t == from_t + g_link_list[link_id].free_flow_link_travel_time)) for from_t in time_set) -        \
                                sum(sum(x[another_agent_id, link_id, from_t, to_t] * from_t for to_t in time_set if (to_t == from_t + g_link_list[link_id].free_flow_link_travel_time)) for from_t in time_set) <= z[agent_id, another_agent_id, link_id] * 10000, "c7")
                        m.addConstr(sum(sum(x[agent_id, link_id, from_t, to_t] * from_t for to_t in time_set if (to_t == from_t + g_link_list[link_id].free_flow_link_travel_time)) for from_t in time_set) -        \
                                sum(sum(x[another_agent_id, link_id, from_t, to_t] * from_t for to_t in time_set if (to_t == from_t + g_link_list[link_id].free_flow_link_travel_time)) for from_t in time_set) >= (z[agent_id, another_agent_id, link_id] - 1) * 10000, "c8")
    print("finish FIFO constraint part 1")
    
    for link_id in link_set:
        for another_link_id in link_set:
            if (another_link_id != link_id) and (g_link_list[link_id].type == 1) and (g_link_list[another_link_id].type == 1):
                for agent_id in agent_set:
                    for another_agent_id in agent_set:
                        if (another_agent_id != agent_id) and (link_id in g_agent_list[agent_id].link_sequence_in_path) and (another_link_id in g_agent_list[another_agent_id].link_sequence_in_path):
                            m.addConstr(z[agent_id, another_agent_id, link_id] == z[agent_id, another_agent_id, another_link_id], "c9")
    print("finish FIFO constraint part 2")   

    

    m.setParam('MIPGap', 0)
    # m.setParam('TimeLimit', 300)
    m.optimize()
    
    
    if m.status == GRB.Status.INFEASIBLE:
        print('Optimization was stopped with status %d' % m.status)
        # do IIS, find infeasible constraints
        m.computeIIS()
        for c in m.getConstrs():
            if c.IISConstr:
                print('%s' % c.constrName)
                    
    # m.write("m.lp")
    
    for agent_id in agent_set:
        from_node_id = g_agent_list[agent_id].from_node_id
        to_node_id = g_agent_list[agent_id].to_node_id
        departure_time = g_agent_list[agent_id].departure_time
        arrival_time = g_agent_list[agent_id].arrival_time
        for from_t in time_set:
            for link_id in link_set:
                for to_t in time_set:
                    current_variable_value = x[agent_id, link_id, from_t, to_t].x
                    if (current_variable_value == 1):
                        # print(x[agent_id, link_id, from_t, to_t, from_s, to_s])
                        g_agent_list[agent_id].link_sequence.append(link_id)
                        g_agent_list[agent_id].time_sequence.append(from_t)
                        g_agent_list[agent_id].node_sequence.append(g_link_list[link_id].from_node_id)
    
    for link_id in link_set:
        from_node_id = g_link_list[link_id].from_node_id
        for from_t in time_set:
            to_t = from_t + g_link_list[link_id].free_flow_link_travel_time
            if (to_t in time_set):
                current_variable_value = y[link_id, from_t].x
                g_link_list[link_id].node_sequence.append(from_node_id)
                g_link_list[link_id].time_sequence.append(from_t)
                g_link_list[link_id].flow_sequence.append(current_variable_value)
                g_link_list[link_id].link_sequence.append(link_id)

    return()                 
                 
                    


def g_write_optimal_solution():
    
    # output passenger file
    workbook = openpyxl.Workbook()
    sheet = workbook.active
    sheet.title = 'output_agent_results'    
    sheet['A1'] = 'agent_id'
    sheet['B1'] = 'node_sequence'
    sheet['C1'] = 'time_sequence'
    sheet['D1'] = 'link_sequence'   
    sheet['E1'] = 'scheduled_departure_time'
    sheet['F1'] = 'scheduled_arrival_time'     
    sheet['G1'] = 'from_node_id' 
    sheet['H1'] = 'to_node_id' 
    sheet['I1'] = 'reservation_flag'
    sheet['J1'] = 'optimized_departure_time'
    sheet['K1'] = 'optimized_arrival_time'
    sheet['L1'] = 'scheduled_total_travel_time'
    sheet['M1'] = 'optimzied_total_travel_time'
    sheet['N1'] = 'arrival_time_deviation'
    sheet['O1'] = 'square_of_arrival_time_deviation'
    
    for agent_id in agent_set:
        row = agent_id + 1
        departure_time = g_agent_list[agent_id].departure_time
        arrival_time = g_agent_list[agent_id].arrival_time
        del(g_agent_list[agent_id].node_sequence[0])
        del(g_agent_list[agent_id].time_sequence[0])
        del(g_agent_list[agent_id].link_sequence[0])
        node_sequence = ";".join(str (node) for node in g_agent_list[agent_id].node_sequence)
        time_sequence = ";".join(str (time) for time in g_agent_list[agent_id].time_sequence)
        link_sequence = ";".join(str (link) for link in g_agent_list[agent_id].link_sequence)
        reservation_flag = g_agent_list[agent_id].reservation_flag
        from_node_id = g_agent_list[agent_id].from_node_id
        to_node_id = g_agent_list[agent_id].to_node_id
        optimized_departure_time = g_agent_list[agent_id].time_sequence[0] - 1
        optimized_arrival_time = g_agent_list[agent_id].time_sequence[-1] + 1
        
        g_agent_list[agent_id].optimized_departure_time = optimized_departure_time
        g_agent_list[agent_id].optimized_arrival_time = optimized_arrival_time
        
        scheduled_total_travel_time = arrival_time - departure_time
        optimzied_total_travel_time = optimized_arrival_time - optimized_departure_time
        arrival_time_devaition = optimized_arrival_time - arrival_time
        square_of_arrival_time_deviation = arrival_time_devaition * arrival_time_devaition
        
        
        sheet.cell(row = row, column = 1, value = agent_id)
        sheet.cell(row = row, column = 2, value = node_sequence)
        sheet.cell(row = row, column = 3, value = time_sequence)
        sheet.cell(row = row, column = 4, value = link_sequence)
        sheet.cell(row = row, column = 5, value = departure_time)
        sheet.cell(row = row, column = 6, value = arrival_time)
        sheet.cell(row = row, column = 7, value = from_node_id) 
        sheet.cell(row = row, column = 8, value = to_node_id) 
        sheet.cell(row = row, column = 9, value = reservation_flag) 
        sheet.cell(row = row, column = 10, value = optimized_departure_time) 
        sheet.cell(row = row, column = 11, value = optimized_arrival_time) 
        sheet.cell(row = row, column = 12, value = scheduled_total_travel_time) 
        sheet.cell(row = row, column = 13, value = optimzied_total_travel_time) 
        sheet.cell(row = row, column = 14, value = arrival_time_devaition)
        sheet.cell(row = row, column = 15, value = square_of_arrival_time_deviation)
        workbook.save('output_agent_results.xlsx')
    
    



    # output link file
    workbook = openpyxl.Workbook()
    sheet = workbook.active
    sheet.title = 'output_link_results'    
    sheet['A1'] = 'link_id'
    sheet['B1'] = 'link_type'
    sheet['C1'] = 'node_sequence'
    sheet['D1'] = 'time_sequence'
    sheet['E1'] = 'flow_sequence'
    row = 1
    for link_id in link_set:
        # if (g_link_list[link_id].type == 1):
        row += 1
        link_type = g_link_list[link_id].type
        node_sequence = ";".join(str (node) for node in g_link_list[link_id].node_sequence)
        time_sequence = ";".join(str (time) for time in g_link_list[link_id].time_sequence)
        flow_sequence = ";".join(str (flow) for flow in g_link_list[link_id].flow_sequence)
        sheet.cell(row = row, column = 1, value = link_id)
        sheet.cell(row = row, column = 2, value = link_type)
        sheet.cell(row = row, column = 3, value = node_sequence)
        sheet.cell(row = row, column = 4, value = time_sequence)
        sheet.cell(row = row, column = 5, value = flow_sequence)

    workbook.save('output_link_results.xlsx')
    
    


    return()




def g_data_statistics():
    
    # calculate cumulative flow counts
    for node_id in node_set:
        g_node_list[node_id].cumulative_departure_count = [0 for t in range(0, g_number_of_time_intervals)]
        g_node_list[node_id].cumulative_arrival_count = [0 for t in range(0, g_number_of_time_intervals)]
        node_type = g_node_list[node_id].type
        if (node_type == 2):
            # origin
            for t in time_set:
                current_cumulative_departure_count = 0
                for agent_id in agent_set:
                    if (g_agent_list[agent_id].optimized_departure_time <= t):
                        current_cumulative_departure_count += 1
                g_node_list[node_id].cumulative_departure_count[t] = current_cumulative_departure_count
        if (node_type == 3):
            # destination
            for t in time_set:
                current_cumulative_arrival_count = 0
                for agent_id in agent_set:
                    if (g_agent_list[agent_id].optimized_arrival_time <= t):
                        current_cumulative_arrival_count += 1
                g_node_list[node_id].cumulative_arrival_count[t] = current_cumulative_arrival_count        

    
    
    
    return()



if __name__=='__main__':
    
    print('Reading data......') 
    g_read_input_data()
    g_generate_in_out_going_link()
    
    start_Gurobi = time.time()
    g_use_Gurobi_to_sovle_the_model()
    end_Gurobi = time.time()
    computing_time_Gurobi = end_Gurobi - start_Gurobi
    print(computing_time_Gurobi)
    
    g_write_optimal_solution()
    g_data_statistics()


    
    